c *******************************************************************
c *                                                                 *
c *        material model # 2 --  deformation plasticity for power- *
c *                               law following initial linear      *
c *                               response with transition region   *
c *                                                                 *
c *                               small-strain analysis only        *
c *                                                                 *
c *******************************************************************
c
c
      subroutine mm02( step, iter, felem, gpn, e, nu, sigyld, exp,
     &                 stress_n, stress_n1, strain, history,
     &                 history1, span, iout, signal_flag )
      implicit integer (a-z)
$add param_def
c
c                   parameter declarations
c   
#dbl      double precision
#sgl      real
     & e(*), nu(*), sigyld(*), exp(*), stress_n(mxvl,*), 
     & stress_n1(mxvl,*), strain(mxvl,*), history(span,*),
     & history1(span,*)
       logical signal_flag
c
c                   locally defined
c
#dbl      double precision
#sgl      real
     & k1, k2(mxvl), l2, pi2, twothd, third, root2, rt23, zero, one,
     & two, half, epsm, e1, e2, e3, epseff(mxvl),
     & epsyld(mxvl), epslim(mxvl), g, c, a, b, 
     & theta, c1, c2, epsnc(mxvl), signc(mxvl), rnc(mxvl),
     & sigeff(mxvl),  eps_elas(mxvl,6), shear_mod(mxvl)
      logical mm02es, debug, signal, nonlinear_flags(mxvl), test,
     &        nonlinear_points
c
      data pi2, twothd, third, root2, rt23
     & / 1.570795, 0.66666667, 0.333333333, 1.4142135623730,
     &   0.816496580927726 /
      data zero, one, two, k1, half
     & / 0.0, 1.0, 2.0, 0.95, 0.5 /  
      data signal / .true. /             
c
c              model history data:
c              (1) 0.0 or 1.0 for linear/yielding
c              (2) scalar effective stress
c              (3) scalar effective strain
c
      debug = .false.
c
      if ( step .eq. 1 ) then
       do i = 1, span
         history(i,1) = zero
         history(i,2) = zero
         history(i,3) = zero
         history(i,4) = zero
         history(i,5) = zero
         history1(i,1) = zero
         history1(i,2) = zero
         history1(i,3) = zero
         history1(i,4) = zero
         history1(i,5) = zero
       end do
      end if
c
c             get the effective strain for the total strain. check if
c             we are still on the linear-reponse below the transition
c             point. set logical linear vs. nonlinear flag for element.
c
      do i = 1, span
         epsm = ( strain(i,1) + strain(i,2) + strain(i,3) ) * third
         e1   = strain(i,1) - epsm
         e2   = strain(i,2) - epsm
         e3   = strain(i,3) - epsm
         epseff(i) = rt23 * sqrt( e1*e1 + e2*e2 + e3*e3 +
     &               half*strain(i,4)**2 + half*strain(i,5)**2 +
     &               half*strain(i,6)**2 )
         epsyld(i) = sigyld(i) / e(i)
         epslim(i) = k1 * epsyld(i) * twothd * (one+nu(i))
         nonlinear_flags(i) =  epseff(i) .gt. epslim(i)
         history1(i,4) = e(i)
         history1(i,5) = sigyld(i)
      end do
c
c             do we have any nonlinear elements at this point 
c             to process.
c
      nonlinear_points = .false.
      nonlin_point     = 0
c
      do i = 1, span
         if ( nonlinear_flags(i) ) then
            nonlinear_points = .true. 
            if ( nonlin_point .eq. 0 ) nonlin_point = i
         end if
      end do
c
      if ( debug ) then
        write(iout,*) ' ' 
        i = max( nonlin_point,1 )
        write(iout,9000) felem+i-1, gpn, e(i), nu(i), sigyld(i),
     &       exp(i), (strain(i,j),j=1,6), 
     &       (stress_n(i,j),j=1,7), history(i,1), history(i,2),
     &       history(i,3), epseff(i), epsyld(i), epslim(i)
      end if
c
c             process points still in linear response region
c
      do i = 1, span
       if ( .not. nonlinear_flags(i) ) then 
        g = half * e(i) / ( one + nu(i) )
        c = e(i) / ( ( one + nu(i) ) * ( one - two * nu(i) ) )
        a = c * ( one - nu(i) )
        b = c * nu(i)
        stress_n1(i,1) = a*strain(i,1) + b*(strain(i,2)+strain(i,3))
        stress_n1(i,2) = b*(strain(i,1)+strain(i,3)) + a*strain(i,2)
        stress_n1(i,3) = b*(strain(i,1)+strain(i,2)) + a*strain(i,3)
        stress_n1(i,4) = g*strain(i,4)                            
        stress_n1(i,5) = g*strain(i,5)                            
        stress_n1(i,6) = g*strain(i,6)  
        stress_n1(i,7) = half * ( stress_n1(i,1)*strain(i,1) +
     &                   stress_n1(i,2)*strain(i,2) +
     &                   stress_n1(i,3)*strain(i,3) +
     &                   stress_n1(i,4)*strain(i,4) +
     &                   stress_n1(i,5)*strain(i,5) +
     &                   stress_n1(i,6)*strain(i,6) )
        history1(i,1) = zero
        history1(i,2) = zero
        history1(i,3) = epseff(i)
       end if
      end do
c      
      if ( .not. nonlinear_points ) go to 1000
c
c
c             process points beyond linear region. iteratively solve for
c             effective stress given effective strain. solution for the
c             effective stress requires a local newtion iteration
c             done by function mm02es. Function is inlined but loop
c             not vectorized.
c
      do i = 1, span
       if ( nonlinear_flags(i) ) then 
        k2(i)  = root2 * ( one-k1 ) / exp(i) + one
        theta  = atan( k2(i)**(one-exp(i))/exp(i) )
        l2     = -tan(pi2-theta)
        c1     = one / (l2 + one)
        c2     = k2(i)**exp(i)
        epsnc(i) = ( two*k1 + l2*c2 - k2(i) ) * c1
        signc(i) = ( two*l2*k1 - l2*c2 + k2(i) ) * c1
        rnc(i)   = root2 * ( k1 - l2*k1 + l2*c2 - k2(i) ) * c1
       end if
      end do
c
      do i = 1, span
       if ( nonlinear_flags(i) ) then 
         test = mm02es( sigyld(i), epsyld(i), nu(i), exp(i),
     &                  epseff(i), sigeff(i), k2(i), epsnc(i),
     &                  signc(i), rnc(i) ) 
         if ( .not. test ) then
           write(iout,9250) felem+i-1, gpn
           call die_abort
           stop
         end if
       end if
      end do
c
c             compute total stresses for total strain. write a
c             yielding message if that option is on.
c             get strain energy  density
c
      if ( signal .and. signal_flag ) then 
       do i = 1, span
        if ( nonlinear_flags(i) .and. history(i,1) .eq. zero ) then 
          write(iout,9300) felem+i-1, gpn, sigeff(i)-sigyld(i)*k1
          history1(i,1) = one
        end if
       end do
      end if
c
      call mm02ns( rnc, signc, k2, k1, exp,
     &             sigyld, epseff, sigeff, nu, e, 
     &             strain, stress_n1, nonlinear_flags, span )
c
c             put effective stress and strain energy into stress
c             vector for point. update history parameters.
c
      do i = 1, span
       if ( nonlinear_flags(i) ) then 
         history1(i,1) = one
         history1(i,2) = sigeff(i)
         history1(i,3) = epseff(i)
       end if
      end do
      if ( debug ) then
        write(iout,*) ' ' 
        i = nonlin_point 
        write(iout,9400) (stress_n1(i,j),j=1,7)
      end if
c
c             compute plastic work density at gauss pts. and
c             the current plastic effective strain.
c         
 1000 continue
      do i = 1, span
         shear_mod(i) = e(i) * half / (one + nu(i))
         eps_elas(i,1) = (stress_n1(i,1) - nu(i)*(stress_n1(i,2)+
     &                   stress_n1(i,3)))/e(i)
         eps_elas(i,2) = (stress_n1(i,2) - nu(i)*(stress_n1(i,1)+
     &                   stress_n1(i,3)))/e(i)
         eps_elas(i,3) = (stress_n1(i,3) - nu(i)*(stress_n1(i,1)+
     &                   stress_n1(i,2)))/e(i)
         eps_elas(i,4) = stress_n1(i,4) / shear_mod(i)         
         eps_elas(i,5) = stress_n1(i,5) / shear_mod(i)         
         eps_elas(i,6) = stress_n1(i,6) / shear_mod(i)         
         stress_n1(i,8) =  stress_n1(i,7) -  half * 
     &    (  eps_elas(i,1) * stress_n1(i,1) 
     &     + eps_elas(i,2) * stress_n1(i,2) 
     &     + eps_elas(i,3) * stress_n1(i,3) 
     &     + eps_elas(i,4) * stress_n1(i,4) 
     &     + eps_elas(i,5) * stress_n1(i,5) 
     &     + eps_elas(i,6) * stress_n1(i,6) )
         stress_n1(i,9) = zero
         if( nonlinear_flags(i) )
     &     stress_n1(i,9) = epsyld(i)*(sigeff(i)/sigyld(i))**exp(i) 
     &                       - sigeff(i) / e(i)
      end do    
c
      return
c
 9000 format('>> debug from mm02. elem, gpn : ',i8,i2,
     & /,    '    e, nu, sigyld, exp : ',4f10.2,
     & /,    '    strains @ n+1:',
     & /,10x,3e15.6,/,10x,3e15.6,
     & /,    '    stresses @ n :',
     & /,10x,3e15.6,/,10x,4e15.6,
     & /,    '    history : ',f4.1,f10.3,e14.4,
     & /,    '    epseff, epsyld, epslim: ',3e14.6 )
 9200 format('    updated elastic stresses: ',
     & /,10x,3e15.6,/,10x,4e15.6 )
 9250 format(//,'>>> Fatal error in mm02',
     &        /,'    mm02es failed to converge for',
     &        /,'    element, point: ',i8,i2,
     &        /,'    job terminated....' )
 9300 format(10x,i8,i3,' point yields.  f = ',f12.3 )
 9400 format('    updated elastic-plastic stresses: ',
     & /,10x,3e15.6,/,10x,4e15.6 )
c
      end
c ********************************************************************
c *                                                                  *
c *   mm02es -- find effective stress for transitional or power-law  *
c *             portions of uniaxial response curves                 *
c *                                                                  *
c ********************************************************************
c
c
      logical function mm02es( sigyld, epsyld, nu, exp, epseff, sigeff,
     &                         k2, epsnc, signc, rnc )
#dbl      implicit double precision ( a-h,o-z)                                
c
#dbl      double precision
#sgl      real
     &   nu, k2
      data  maxitr / 30 /, toler / 0.0001 /
      data  twothd, third, one, two, half, zero 
     &  / 0.666666667, 0.333333333, 1.0, 2.0, 0.5, 0.0 /
c
c             use newton-raphson iterative solution to find the
c             effective stress given the effective strain for
c             a multi-dimensional deformation model. the uniaxial
c             curve is a linear region, a circular transition
c             region and a power-law region.
c
c             sigyld  -- yield stress
c             epsyld  -- yield strain corresponding to sigyld
c             nu      -- poisson's ratio
c             exp     -- power-law exponent
c             epseff  -- effective uniaxial strain (given)
c             sigeff  -- effective uniaxial stress (computed)
c             mm02es  -- .true. if procedure converges
c             k2      -- upper point on transition
c             epsnc
c             signc   -- center of circle for transition (given)
c             rnc     -- radius of transition circle (given)
c
c
      iterno = 1
      mm02es = .false.
      expm1  = exp - one
      trnlim = ( k2**expm1 + twothd * ( nu - half ) ) * k2 * epsyld
      epsn   = epseff / epsyld     
c
c             handle circular transition and power-law in
c             separate iteration loops
c
      if ( epseff .ge. trnlim ) go to 200
c
c             perform newton raphson on circular transition part of
c             response to compute effective stress. starting guess for
c             normalized effective stress is k2 point on transition curve.
c             for very low hardening, e.g. n=20, with an effective strain
c             very close to the upper transition, the iteration process
c             will blow-up. this is detected by checking the sign of
c             radical. if a blow-up is detected, simply switch to the power-law
c             part of curve and convergence will be quick!
c              
      const = ( two * nu - one ) * third
      sign  = k2
      radical = sqrt( rnc*rnc - (sign - signc)**2 )
      resid   = epsn - epsnc + radical - const * sign
 100  continue                            
      dfdsn = (signc - sign ) / radical - const
      signew = sign - resid / dfdsn
      radical = rnc*rnc - (signew - signc)**2 
      if ( radical .lt. zero ) then
         iterno = 1
         go to 200
      end if
      radical = sqrt( radical )
      resid   = epsn - epsnc + radical - const * signew
      if ( abs(signew-sign) .le. toler ) then
        sigeff = signew * sigyld
      else
        sign = signew
        iterno = iterno + 1
        if ( iterno .gt. maxitr ) return
        go to 100
      end if
      mm02es = .true.
      return
c
c             perform newton raphson on power-law part of
c             response to compute effective stress. starting guess for
c             normalized effective stress is n th root of normailzed
c             effective strain (same approach used for regular deformation
c             plasticity model).
c
 200  continue                           
      sign  = epsn**(one/exp)
      c     = twothd * ( nu - half )
      resid = sign**exp + c*sign - epsn
 210  continue
      dfdsn  = exp*sign**expm1 + c
      signew = sign - resid / dfdsn
      resid  = signew**exp + c*signew - epsn
      if ( abs(signew-sign) .le. toler ) then
        sigeff = signew * sigyld
      else
        sign = signew
        iterno = iterno + 1
        if ( iterno .gt. maxitr ) return
        go to 210
      end if
      mm02es = .true.
      return
c
      end
c ********************************************************************
c *                                                                  *
c *   mm02ns -- new stresses and strain energy density               *
c *                                                                  *
c ********************************************************************
c
c
      subroutine mm02ns( rnc, signc, k2, k1, exp,
     &                   sigyld, epseff, sigeff, nu, e,
     &                   strain, sig, nonlinear_flags, span )
      implicit integer (a-z)
$add param_def
c
c             parameter declarations
c
#dbl      double precision
#sgl      real
     & rnc(*), signc(*), k2(*), k1, exp(*), sigyld(*),
     & epseff(*), sigeff(*), nu(*), e(*), strain(mxvl,*), sig(mxvl,*)
      logical  nonlinear_flags(*)
c
c             local declarations
c
#dbl      double precision
#sgl      real
     & one, two, three, half, onefive, six, pt25, twothd, 
     & third, const1, const2, epsv, sigm, upress(mxvl),
     & sigk1, uelask1(mxvl), sncso(mxvl), rncso(mxvl),
     & c1(mxvl), rncsq(mxvl), betak1, uk1(mxvl), term1,
     & term2, term3, betase, uk2, sigk2, betak2, expp1, c2,
     & term4, term5, term6, term7 
c
      data one, two, three, half, onefive, six, pt25
     &   / 1.0, 2.0, 3.0, 0.5, 1.5, 6.0, 0.25 / 
      data twothd, third
     &  / 0.6666666667, 0.33333333333 /                         
c
c             compute total stress given total strain and
c             effective strain. also the total work density.
c
c
c             compute stress by volumetric and deviatoric
c             components.
c
      do i = 1, span
       if ( nonlinear_flags(i) ) then 
        const1 = twothd * sigeff(i) / epseff(i)
        const2 = e(i) / ( one - two*nu(i) ) - const1
        epsv   = ( strain(i,1) + strain(i,2) + strain(i,3) ) * third
        sig(i,1) = const1 * strain(i,1) + const2 * epsv
        sig(i,2) = const1 * strain(i,2) + const2 * epsv
        sig(i,3) = const1 * strain(i,3) + const2 * epsv
        sig(i,4) = const1 * strain(i,4) * half
        sig(i,5) = const1 * strain(i,5) * half
        sig(i,6) = const1 * strain(i,6) * half
       end if
      end do
c
c             compute exact strain energy density. we have a point on
c             transition part of curve or power-law part of curve.
c             energy up to k1*sigyld is same for both cases.
c
      do i = 1, span
       if ( nonlinear_flags(i) ) then 
        sigm      = third * ( sig(i,1) + sig(i,2) + sig(i,3) )
        upress(i) = ( onefive - three*nu(i)) * sigm * sigm   
        sigk1     = k1 * sigyld(i)
        uelask1(i)= third * (one+nu(i)) * sigk1 * sigk1
        sncso(i)  = signc(i) * sigyld(i)
        rncso(i)  = rnc(i) * sigyld(i)
        c1(i)     = ( two*nu(i) -one ) / six
        rncsq(i)  = rncso(i) * rncso(i)
        betak1    = asin( (sigk1 - sncso(i))/rncso(i) )
        term1     = c1(i) * sigk1 * sigk1
        term2     = -sigk1 * sqrt( rncsq(i) - (sigk1-sncso(i))**2 )
        term3     = rncsq(i) * ( half*betak1 + pt25*sin(two*betak1) )
        uk1(i) = term1 + term2 + term3
       end if
      end do
c
      do i = 1, span
       if ( nonlinear_flags(i) ) then 
        if ( sigeff(i)/sigyld(i) .lt. k2(i) ) then
c
c                 Stress point in transition:
c                   1) pressure component (total)
c                   2) deviatoric on linear curve up to k1*sigyld
c                   3) deviatoric transition from k1*sigyld to sigeff.
c
        betase   = asin( (sigeff(i) - sncso(i))/rncso(i) )
        term1    = c1(i) * sigeff(i) * sigeff(i)
        term2    = -sigeff(i) * sqrt( rncsq(i) -
     &              (sigeff(i)-sncso(i))**2 )
        term3    = rncsq(i) * ( half*betase + pt25*sin(two*betase) )
        uk2      = term1 + term2 + term3
        sig(i,7) = ( upress(i) + uelask1(i) + (uk2-uk1(i)) ) / e(i)
c
      else
c                 Stress point in power-law region:
c                   1) pressure component (total)
c                   2) deviatoric on linear curve up to k1*sigyld
c                   3) deviatoric transition from k1*sigyld to k2*sigyld
c                   4) deviatoric power-law from k2*sigyld to sigeff
c
        sigk2    = k2(i) * sigyld(i)
        betak2   = asin( (sigk2 - sncso(i))/rncso(i) )
        term1    = c1(i) * sigk2 * sigk2
        term2    = -sigk2 * sqrt( rncsq(i) - (sigk2-sncso(i))**2 )
        term3    = rncsq(i) * ( half*betak2 + pt25*sin(two*betak2) )
        uk2      = term1 + term2 + term3
c
        expp1    = exp(i) + one
        c2       = (one/sigyld(i)**(exp(i)-one)) * exp(i) / expp1
        term4    = c1(i) * sigk2 * sigk2
        term5    = c2 * sigk2**expp1
        term6    = c1(i) * sigeff(i) * sigeff(i)
        term7    = c2 * sigeff(i)**expp1
        sig(i,7) = ( upress(i) + uelask1(i) + (uk2-uk1(i)) +
     &              (term6+term7) - (term5+term4) ) / e(i)
c
        end if
       end if
      end do
c
      return
c
      end

c     ****************************************************************
c     *                                                              *
c     *                 subroutine mm02_set_sizes                    *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified: 12/14/14 rhd                *
c     *                                                              *
c     *    called by warp3d for each material model to obtain        *
c     *    various sizes of data for the model                       *
c     *                                                              *
c     ****************************************************************
c
      subroutine mm02_set_sizes( info_vector )
      dimension info_vector(*)
c
c        set infor_data
c
c         1        number of history values per integration 
c                  point. Abaqus calles these "statev". Values
c                  double or single precsion based on hardware.
c    
c         2        number of values in the symmetric part of the 
c                  [D] for each integration point. for solid
c                  elements this is 21, for cohesive elements this 6.
c
c         3        = 0, the material model returns "unrotated"
c                       Cauchy stresses at n+1
c                  = 1, the material model returns the standard
c                       Cauchy stresses at n+1
c
c         4        number of state variables per point to be output
c                  when user requests this type of results
c
      info_vector(1) = 5
      info_vector(2) = 21
      info_vector(3) = 0
      info_vector(4) = 3
c
      return
      end
c
c     ****************************************************************
c     *                                                              *
c     *             subroutine mm02_states_values                    *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *               last modified : 1/9/15 (rhd)                   *
c     *                                                              *
c     ****************************************************************
c
      subroutine mm02_states_values( itype, elem_states_output,
     &                                nrow_states, num_states  )
c
c                       access some global data structures
c
      use elem_block_data, only: history_blocks, history_blk_list
      use main_data, only: elems_to_blocks
c      
      implicit integer (a-z)
$add common.main
c
c                       parameters
c
      integer :: nrow_states, itype, num_states
#dbl      double precision :: elem_states_output(nrow_states,*)
#sgl      real  :: elem_states_output(nrow_states,*)
c
c                       locals
c
#dbl      double precision, 
#sgl      real,
     & allocatable :: history_dump(:,:,:), one_elem_states(:)
      integer :: relem, elnum, hist_size, blockno
      logical :: do_a_block
      double precision :: zero
      data zero / 0.0d00 /
c      
c           build deformation plasticity states values output.
c
c              itype > 0 => this is the block number. do all elements
c                           in the block
c
c              itype < 0 => this is an element number. put state
c                           values into column 1 of results.
c 
      do_block = .true.
      if( itype. gt. 0 ) then
         do_a_block = .true.
         blockno = itype
      else
         do_a_block = .false.
         elnum = -itype
         blockno = elems_to_blocks(elnum,1)
      end if          
c
      local_debug = .false.
      felem       = elblks(1,blockno)
      elem_type   = iprops(1,felem)
      mat_type    = iprops(25,felem)
      int_points  = iprops(6,felem)
      span        = elblks(0,blockno)
      hist_size   = history_blk_list(blockno)
      
      if( local_debug ) write(out,9050) blockno, felem, elem_type,         
     &         mat_type, int_points, span, hist_size
c
c           temporary block of history so it can be re-organized
c
      allocate( one_elem_states(nrow_states) )
      allocate( history_dump(hist_size,int_points,span) )
      history_dump = reshape( history_blocks(blockno)%ptr,
     &           (/hist_size,int_points,span/) )
c      
      if( do_a_block ) then    
        do relem = 1, span
           elnum = felem + relem - 1  ! absolute element number
           one_elem_states(1:nrow_states) = zero 
           call mm02_states_values_a
           elem_states_output(1:nrow_states,relem) = 
     &                one_elem_states(1:nrow_states)
        end do
      else
        relem = elnum + 1 - felem
        one_elem_states(1:nrow_states) = zero
        call mm02_states_values_a
        elem_states_output(1:nrow_states,1) =
     &                one_elem_states(1:nrow_states)
      end if  
c        
      deallocate( history_dump, one_elem_states )
c        
      return
c      
 9050 format(10x,"block, felem, etype, mtype:  ",4i7,
     &  /,10x,   "int_pts, span, hist_size:    ",3i7 )
c 
      contains
c     ========      
     
c     ****************************************************************
c     *                                                              *
c     *                 subroutine mm02_states_values_a              *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *                   last modified : 1/9/15 (rhd)               *
c     *                                                              *
c     ****************************************************************
c
      subroutine mm02_states_values_a
c   
      implicit none      
c
c                       locals
c
      integer :: ipt   
#dbl      double precision :: 
#sgl      real :
     & nl_flag, sig_eff, eps_eff, ymod, sigyld, epsyld
c
      nl_flag = zero
      sig_eff = zero
      eps_eff = zero
c       
      do ipt = 1, int_points
        ymod    = history_dump(4,ipt,relem)
        sigyld  = history_dump(5,ipt,relem)
        epsyld  = sigyld / ymod
        nl_flag = nl_flag + history_dump(1,ipt,relem)
        sig_eff = sig_eff + history_dump(2,ipt,relem) / sigyld
        eps_eff = eps_eff + history_dump(3,ipt,relem) / epsyld
        
      end do
c
      one_elem_states(1) = nl_flag / dble(int_points)            
      one_elem_states(2) = sig_eff / dble(int_points)
      one_elem_states(3) = eps_eff / dble(int_points)            
c
      return     
c
      end subroutine mm02_states_values_a
      end subroutine mm02_states_values
c     ****************************************************************
c     *                                                              *
c     *                 subroutine mm02_states_labels                *
c     *                                                              *
c     *                       written by : rhd                       *
c     *                                                              *
c     *               last modified : 1/11/2015 (rhd)                *
c     *                                                              *
c     ****************************************************************
c
      subroutine mm02_states_labels( size_state,
     &      num_states, state_labels, state_descriptors, out,
     &      comment_lines, max_comment_lines, num_comment_lines )
      implicit none
c
c                       parameters
c
      integer :: size_state, num_states, out, max_comment_lines,
     &           num_comment_lines
      character(len=8)  :: state_labels(size_state)
      character(len=60) :: state_descriptors(size_state)
      character(len=80) :: comment_lines(max_comment_lines)
c
c                       locals
c
      integer :: i
      logical, save :: do_print = .false.
c
      num_states = 3     
      state_labels(1) = "nl flag"
      state_labels(2) = "seff/s0"
      state_labels(3) = "eff/e0"
      state_descriptors(1) = "=0 all pts linear, =1 all pts nonlinear"
      state_descriptors(2) = "mises / sig_yld"
      state_descriptors(3) = "eps-eff / eps_yld"
c
      num_comment_lines = 0
c      
      if( do_print ) then
        do i = 1, 3 
          write(out,9010) i, state_labels(i), state_descriptors(i)
        end do
        do_print = .false.   
      end if
c          
      return
 9010 format(2x,i3,2x,a8,2x,a)      
      end